o 

(N 



> 



Perfect simulation for the Feynman-Kac law on the path 

space 

Christophe Andrieu* Nicolas Chopin f Arnaud Doucetf Sylvain Rubenthaler § 

2nd October 2012 



O 

Abstract 

I This paper describes an algorithm of interest. This is a preliminary version and we intend 

on writing a better descripition of it and getting bounds for its complexity. 

5? 

1 Introduction 

We are given a transition kernel M (on a space E) , Mi a probability measure on E and potentials 
{Gk)k>i (Gk ■ E — » M+). We want to draw samples according to the law (on paths of length P) 

E(f(X 1 ,...,X P )il^ 1 1 GdX,) 
f — , where (Xfe) is Markov with initial law Mi and transition M. For all n <G N, we note [n] = {1, . . . , n}. 

m 
o 

2 Densities of branching processes 
CN 2.1 Description of a branching system 



We start with TVi particles (i.i.d. with law Mi, iVi is a fixed number). If we have iV* particules at 
time n, the system evolves in the following manner: 

• The number of childern of X l n (the i-th particle at time n) is a random variable A l n+1 with 
law f n +i such that : P(yl5j +1 = j) = f n+ i(G n (X^),j) (here, /„ is a law with a parameter 
G n (X l n ), we will define this law later). The variables A l n+1 (1 < i < N n ) are independent. 
We than have N n+1 = 4 +1 

• We draw a n+ i uniformly in S]y n+1 (the jV n+ i-th symmetric group). 

. We set Vj G [iV„], B j n+1 = {A l n+1 + • • • + A{- + \, . . . , ^ +1 + • • • + ^ + If i e 

we draw ~ M(X^ .). 

Such a system has a density on the space 
{(n 2 , . . . ,np,x l n ,A l n ,<T n ) : n 2 , ■ ■ ■ ,n P e N, 
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x i n £ E(l<n<P,l<i< n n ),A i n eN(2<«<F,l<!< n n ),a n £ S N J2 < n < P)} . 
This density is equal to : 

qo(N 2 , ■ ■ ■ , Np, (A l n )2<n<P,l<i<N n , (Xn)l< n <P,l<i<N n , {&n)2<n<p) 

JVi P JV n _i 

»=1 n=2 i=l jeT n {B' n ) 

The random permutations ctat ease the writing of the formulas but have no deep signification. 

2.2 Proposal density 

We take the above branching system and we draw a path by drawing a number i uniformly in 
{1, . . . , N P } and taking the path of the ancestors of X l p . The branching system plus this trajectory 
live in the following space 

{(n 2 , . . . , n P , x l n , A % n , a n , h) : n 2 , . . . , n P £ N, 

<e£(l<n<P,l<i< n n ),^ e N(2 < n < P, 1 < i < m), 

<J n eS Nn (2<n<P),b t e [Ni](l < i < n)} , (2.1) 

and have the following density : 

q(N 2 , ■ ■ ■ ,Np, (A l n )2<n<P,l<i<N n , i x n)l<n<P,l<i<N n , {&n)2<n<Pi (bk)l<k<p) 

= T^9o(^2, • • • , N P , (A I n ) 2 < n <P,l<i<N n , (Xn)l< n <P,l<i<N n , {^n)2<n<p) ■ 
iv p 

2.3 Target law 

We draw a trajectory (j/i,. . . ,yp) with the law 7r then a branching system conditioned on con- 
taining the trajectory (yi, . . . , j/p). The order of operations is as followed 

• Draw (j/i, . . . , yp) with law tt(.). 

• We draw £>i uniformly in [Ni], we set x^ 1 = y 1 . We draw (£i)i<i<jVi,i^6i i-i-d. variables of 
law Mi. 

• If we have the (n — l)-th generation, we draw A^T^ 1 with law /(G„_i(x^J 1 1 ), .) conditioned 
to be in N* (we call this law f(G n -i(x r ^Si), .)). For i £ iV„_i, i 7^ we draw <~ 
/n(G„_i(a;^ 1 1 ), .). We set iV n = X^Sf 1 A l n . Weaw a n uniformly in <Sjv„- We set b n = 
o n (A\ + ■■■ + A'n- 1 - 1 + 1), afr = j/„. For j e [JV n ], if j ^ b n and j G u„(S^) (B* - 

+ • • • + ^T 1 + 1, •■•,< + ••• + 4,}), we draw < ~ M«_ 1; .). 

We get a variable in the following space 

{(n 2 , . . . ,n P ,x l n ,A z n ,a n ,bi) : n 2 , . . . ,n P £ N*,x % n £ E(l < n < P,l < i < n n ), 

A\ £ N(l < n < P,l < i < n n ),a n £ S N J2 <n< P),bi £ [Ni](l < i < n)} , 

with the following density: 

9(N 2 , ■ ■ ■ , Np, (A l n ) 2 < n <PA<i<N n , (Xn)l< n <P,l<i<N n , (o"n) 2 <n<P, (&fe)l<fe<p)) 

=^,...,4-) A- n 

1 l<i<N 1 ,i^b 1 
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n=2 V l<i<Af„_i,i^&„_i 

x ^t n n ■ ^ 

i<j<iv„_i je<r n (Bi)„jjtb n J 
Notice that: (Vz, fc) f n {g,k) = ( £ n 1i is conditioned on having at least one children). 

2.4 Ratio of the densities 

We write the ratio n/q and we get: 

n(N 2 , ■ ■ ■ , N P , (A£ [ )i< n <p_i j i<j<jv n , (a4)l<ri<P,l<i<AT„ , (<J n )2<n<P, (bk)i<k<p)) 
q(N 2 , ■ ■ ■ , Np 7 (4jj)l< n <p-l,l<i<iV n , (^?j)l<n<P,l<i<iV„, (&n)2<n<P, (bk)l<k<p)) 

, bl bp s v 1 /n^n-l^n-l),^"- 1 ) 

= 7TIX! 1 , ...,I„M X — — X 5 r r X r r . 

Nl M 1 (x^)Un=2M(x b n ^^ h n) n=2 /» (G„-l ) , 1 ) 

Let us take /„ such that for all (i ^ 0), y^fe] = ^ f° r some comstant /?„. This means that 
l-/n(<7,0)= We then get: 

,p 



7T 



(•••)_ ^Fnr= 2 A 



<?(...) 

with Z = E(n„=i G n {X n )) ((X n ) „>X is a Markov chain with initial law M\ and kernel transition 
M). 



3 Perfect simulation algorithm 

3.1 Stability of the branching process 

We want the branchin process to be stable. So we need that 

JV„_i +oo 

]vT" ^2 E^^^i)-^ be of order 1 ( Vn )- ( 3J ) 

™ 8=1 j = l 

Let us take: /3„ > 1 1 1 1 e» (Vn), and (for some fc n ), f n (g,0) = 1 - -jfc, f n {g,i) = 7^ pour 
1 < i < fc„. We then get Yli=i * x /n(.9i *) — gg+llg _ g G ft j s sensible to fix fc„ such that 

k 4-1 1 W 

^» = ^ ti x^G„_ 1 «_ 1 ) (3.2) 

i=l 

where is a sequential Monte-Carlo system with N particles, this has to be computed 

beforehand. Simulations show that this procedure indeed gives you stable branching processes. 

3.2 Markovian transition 

We know want to use a backward coupling algorithm (as in [FT98 PW96J). The integer Ni is 
fixed. We take (z\, . . . , Zp) E E p . 
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• We draw iV2, . . . ,Np, (^^Kn^RKK^.^Bn, (j4^)l<n<P,l<j<n„> (£>n G <SjV„ )2<n<P , (Pfc)l<fc<P 

with the density 



7r(z 1; . . . ,z P ) 



(zi,...z P in place of equation fl2.2|). This amounts to drawing a genealogy 



conditionned to contain (z\, . . . , Zp). Let us set Vn G {1, . . . , P}, ™ = z n . Let A" be the 
variable containing all the N n , X l nl A l n , S n , B n . 

• WedrawiV 2 , ...,N P {X^^^p^^, {A l n )i< n <p,i<i<N^ (S n & S Nn ) 2 <n<p, [B k )i<k<p 
with density q(.). We denote by X the corresponding variable. 

• With probability inf (l, §p^7§} ) , we set (Z u . . . , Z P ) = (Zf \ . ..,X P Bp ), and with the 
complementary probability, we set (Z±, . . . , Zp) = {z\, ■ ■ ■ , Zp). 

The transformation of (zi, . . . , zp) into (Zi, . . . , Zp) is a Metropolis Markov kernel (on E p ) for 
which 7r is invariant (much in the spirit of [ADH10 ). Recall that 



v(X)q(X) Np 



(3.4) 



3.3 Backward coupling 



We are given i.i.d. variables (Uq, £/_i, U-2, ■ ■ ■ )■ Any U-i is sufficient to make a simulation 
of the Markovian transition above. We introduce a function F parametrizing this transition 
(we can write the transition in the following manner: (Z±, . . . , Z P ) = Fxj{z\, . . . , zp)). By 
Theorem 3.1 of [FT98 . if T is a stopping time, relatively to the filtration [ct{Uq, . . . , J7_i)),>o, 
such that Vfzj 1 ',. ■ • , 4 1} ). (4 2 \ 4 2) ) e E p , F V _ T o • • • o ^(z^, . . . , z^) = F, 
Fu {zf \ . . . , Zp" 1 ), then F[/_ T o • • • o FyJ zj, . . . , Zp^) is exactly of law tt. 
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We now look for a lower bound of (3.4 1 for a trajectory [z%, . . . , Zp) G E and i G N, ?7_j fixed. 
We add the following hypothesis. 

Hypothesis 1. There exists a function / : K — > R such that: for all x\ G E, i G N, fixed, 
(xi, . . . ,xp) trajectory drawn with transitions M using the variables [/_; (which we will denote 
by x j+1 = Mc/_ j (x,), Vj G {2, . . . , P}), for all Se > 0, Vj G {2, . . . , P}, diam(M° (j - 1) (P £ (xx))) < 
/°a-D(e). 

Example 3.1. If the transition M is (for some constants a,b) : 

1 / (w 

M(x, dy) = —j= exp ' 



V2^ * V 2^ 
then we can take fix) = ax. 

We now want to bound the number of descendants generated by the trajectory (zi,... ,Zp) 
during the conditional drawing using the variables U-i. Let us precise how we do this conditional 



drawing (zi, . . . , 2p). We fix Vn, /3 n = ||(rn||oo and k n satisfying (3.2 K For g G [0;||G||oo], we 
set u i— > F~g(u) to be pseudo-inverse of the cumulative distribution function of the law f n (g, .) 
and we set u i— > F~g(u) to be the pseudo-inverse of the cumulative distribution function of the 

law f n (g,.). We are given a family (V u , W u ) ue (m*)["] «>i (random variables indexed by infinite 
sequences of N*) of independent variables of law U([0; 1]). We are given (cr n ,N)n,N>i independent 
variables such that \/n,N. a n .N is uniform in Sn- Suppose there exists M' : E x [0; 1] —> E such 
that if U ~ U{[0; 1]), x G E then M'(x, U) ~ M(x, dy). Suppose there exists M{ : [0; 1] — >• M tsuch 
that if U ~ W([0; 1]), then M((J7) - Mi(dx). The simulation goes as follows. 

• We set Xf = A/{ (Vm) ((i) is a sequence of length 1 taking value i) for all i G [A r i]\{l}, and 
X\ =zi. We define *i : [JVi] (N*) [1] by * x (i) = (i). 
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Suppose we have made the simulation up to time n < P and we have a function : [N„] — > 
(N*)^"! (describing the genealogy of the particles. & n (i) is the complete ancestral line of 
particle i). 

- For i G [N n ]\{l}, we set A\ +1 = F', l (x ^(W 9n(i) ); 



and if i = 1 , then X\ = z n , 



and we set A l n+1 = ^ n + 1)Gn(zn) (W* n(i) ). We set N n+1 = A^ +1 . 

- For j G [JV B+1 ]\{1}, if A\ +1 + ■■■ + A l - + \ < j < A\ +1 + ■■■ + A l n+Il we set * n+1 (j) = 

- and if j = 1. we set = z n+1 , $ n +i(j) = (1, 1, . . . , 1). 

We then set X[ = x^ Nl(i) {\<i< N x ), B x = cr^(l) (beware, B; and Bj have different 

meanings). We then proceed by recurrence. If we have (^j)i<j<n,i<i<N n -, {Aj)2<j<n,i<i<Nj-n 

(o-j)2<i<n, B 1 ,...,B n with X;. - xJ 3% ^ W (Vj G [n], i G [ATj]) then: 

We set A n+1 = A n ^f n , B' n+1 = + • ■ ■ + + 1, . . . , + ■ • • + 

<T n +i = K+i = o-n+i(^n+f nW ): = JC£ 1 '*' B+1 * , (Vi. . .)• We have 

- if i G B* = <r~l 1>N (K+f " ) and * ^ B„+i := ^nji.jv^C 1 )) °"n+i,iv„ +1 («) e 

S n+f* * > X n+l' N " +1 = Af '(- X 'n n ' JVn(? >^*»(<r„,jv n («)))) tnen X n+1 = M ' ' [ X ni V 9 n {<T n ,N n («))) 



— and in the case i = B n+ i, A^"^ 1 = X^+i = z n +i 
And we have 



- if B n+1 i B l n+1 , then #B^ +1 = #fl£f" W = <-f» W = W^,^,)), 

- if B n+1 G St +lj then - /^M^,^,))- 



This procedure draw (X nl A n , B n , a n ) with the density (3.3 1 (in pratice, one can get rid of the 
simulation of the permutations since they have no influence on the trajectories we are interested 
in). We will note (X l n , A l n , B n ,<r n ,n £,..) = $((z 4 ) 4e[P] , (V u , W u ) ue(N , )H ,„>i , (G n )i< n <p)- 

Lemma 3.2. If in the procedure above, we replace A l n+1 = F~ +1 G , z \{W^ n {i)) (in the case 

*n(i) = (Ni, 1, • ■ • , 1)) by A l n+1 = F~l 1H ^ Zn) (W^, n (i)) for some function H„ > G n , then we get 
a different system, which we will note with~ 

(X;,I;,S„,2F„,n G . . . ) = *((3{) ig [j>], (K, Wu) u6(N . )Win >i , (H n )i< n <p)) , 

such that Vn, 1 < n, < iV„} C {X l n , 1 < n < A^}. moreover, the descendants of Z\, . . . ,Zp at 
time P are independent variables. 

Let 6 > 0. For all n G [P], let us take Pi = G±, . . . , H n -i — G n -i t and for k > n. 

Hk{x) = (sup ly _ Zkl<f ^- nHs) G n (y) if \x - 2*| < /°( fe -")(J) 
1 G„ (y) otherwise , 

and let us note with ~ the corresponding system, 

meaning (X* n , A\,B n ,a n ,n€...) = $((zi) ie[P] , (V u , W u ) ue(N «)H in >i , (P n )i<n<p)) ■ 
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Let z[, . . . , z' P be such that z[ € Bg(zi), Vi. We have 

(A^,A l „,5„,a„,ne ...) = $((^ e[ P],(U u ,IU u ) ue 

(N*)H,n>l ! (Cn)l<n<p) • 

Using the Lemma above and Hypothesis [l] we have Np < Np. Let $' be such that 
N P = &((zi) ie[P] ,(V u ,W u ) ue 

(N*)W,n>l ' {H n )l<n<p)) ■ 

3.4 Examples 

3.4.1 Gaussian example 

We draw sequences (X n ) ne [ P ], (Y n ) ne[P] such that: X\ ~ 7V(0, 1), X n+1 = aX n + bV n+1 (a e]0; 1[), 
y„ = X„ + cW n with i.i.d. variables U„, W n of law 7V(0, 1). We set 

Mi(dx) = -J=e~ x2 / 2 dx, M(x,dy) — exp ^— ^ff^ J dy. We want to bound, at time P, 

the particles descending from a fixed trajectory. The descendants of different z n are independant 
so we look, for all n, at which is the z n spawning the most descendants at time P. Using the 
result above, we slice E in balls of size S. If z' n is in a ball of size S containing z n , the number of 
descendants of z' n at time P computed with potentials G. is bounded by the number of descendants 
of z n at time P computed with potentials H . The potentials G n going to at ±cx), we do not 
have to explore the whole of M, as soon as z n is far enough from Y n so that it has children under 
potential H n , we can stop the exploration. 

Remark 3.1. With 5 = (or S very small), if we look at the number of descendants at time P 
of an individual at time n and we maximise in the position of the individual, we will finite some 
finite quantity (not exploding when P — n — > +oo. For the maximisation step, we have to take 
6 > and then this maximum explodes (slowly). So, there a balance to find between 5 small 
(maximisation step takes a lot of time) and 5 big (explosion in the number of particles) . A rule of 
thumb, coming from the experience, is that the population do not explode as long as the number 
of children per individual is of order 2, 3. 

3.4.2 Directed polymers 

Let (X n ) n >i be a symmetric simple random walk in Z with X\ = 0. We are given i.i.d. variables 
(£n,i)n>i,iez with Bernoulli law of parameter p > 0. We set (random) potentials : V n (i) = 
exp(— /3£„ (/3 > 0) and we are interested in the following law (quenched, meaning the £„ j are 
fixed) : 

E C (/(x 1:n )nLi^(^)) 

TlmtfJ- E^UUV k (X k )) • 

This kind of model is described in |BTV08| . If we take the expectation over all the variable: 
E(max de la traj. sous 7Ti :n ) behaves as rfi with ( ^ 1/2. 

Using our algorithm, we can simulate trajectories under the law nup (for fixed £, P € N*). 
The research of the ancestors having the biggest number of descendants at time P makes that the 
computational cost is P 2 . Here is the drawing of E(max . . . ) as a function of n in a log- log scale 
(the blue line has gradient 2/3, the green line has gradient 1/2): 
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Figure 3.1: gradient estimation (least square)=0,63 
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